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C""^ ■ So far the problem of a spin glass on a Bethe lattice has been solved only at the 

^^ ' replica symmetric level, which is wrong in the spin glass phase. Because of some technical 

difficulties, attempts at deriving a replica symmetry breaking solution have been confined 

(~| ' to some perturbative regimes, high connectivity lattices or temperature close to the critical 

f~| I temperature. Using the cavity method, we propose a general non perturbative solution 

of the Bethe lattice spin glass problem at a level of approximation which is equivalent to 

a one step replica symmetry breaking solution. The results compare well with numerical 

simulations. The method can be used for many finite connectivity problems appearing in 

combinatorial optimization. 
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I. INTRODUCTION 



T3 

O ' The spin glass problem has been around for twenty five years, but its understanding has turned out 

, ^, to be remarkably complicated. It is generally considered as solved only in its fully connected version 

introduced by Sherrington and Kirkpatrick [|l| . The first consistent solution was derived with the replica 

method |2,H] and it was then confirmed using a probabilistic approach, the cavity method, which avoids 

K*" ' the strange (and powerful) mathematical subtleties of the replica approach y,|[ . A rigorous proof of the 

00 ' validity of the solution is still lacking, in spite of recent progress P-R]. 

A slightly more realistic theory of spin glasses, still of the mean field type, deals with the situation in 

which each spin interacts only with a finite number of neighbours. Models of this type include the spin 

^^' t glass on a Cayley tree, a Bethe lattice and a disordered random lattice with fixed or with fluctuating 

f~^ connectivity. There are many motivations for studying such problems. On one hand one may hope 

f^ J to get a better knowledge on the finite dimensional problem, since these models include a notion of 

neighborhood which is absent in the infinite range case. But another motivation is the possibility to solve 

C^ these problems using different methods, like iterative methods which are typical of statistical mechanics 

on tree-like structures. In fact the cavity method is a generalization of the Bethe Peierls iterative method 

to the case in which there may exist several pure states, and it is therefore very natural to work out 

the details of this generalisation, and to test its validity. Another important aspect comes from the 

connection between the statistical mechanics of disordered systems and the optimization problems: many 

^ , of the interesting random optimization problems turn out to have a finite connectivity structure. This is 

the case for instance of the travelling salesman problem Isl , the matching the graph partitioning |10| , pl| 

or the K-satisfiability problem Jl^ l. 

While the problems of a spin-glass on tree-like lattices were naturally studied very soon after the 
5-H ' discovery of spin glasses, the present status of the knowledge on these systems is still rather poor compared 
to that on the SK model. A lot of efforts have been devoted to the simple Bethe Peierls method which 
builds up a solution in terms of the distribution of local magnetic fields [0-18|. However this simple 
iterative solution, which may be relevant for a Cayley tree with a certain type of boundary conditions 
p8[ is wrong for the Bethe lattice spin-glass. When the replica formalism is used, this simple iterative 
solution turns out to be equivalent to the replica symmetric (RS) solution. However one knows that there 
exists a replica symmetry breaking (RSB) instability similar to the one found in the SK model |l^,n9-p3]. 
Unfortunately, in the replica formalism, the RSB solution could be found only in some rather limited 
regimes: expansion around the high connectivity (SK-like) limit p3], or close to the critical temperature 
p(|. The main problem encountered in all these attempts is a very general one, common to all disordered 



problems with a finite connectivity. Roughly speaking it can be summarized as follows: the distribution 
of local fields, even within one pure state, is not a simple Gaussian as in the infinite range problems, but 
a more complicated function [g,|4 25 1. When one takes into account the existence of several pure states. 



the natural order parameter, even within the simplest "one step" replica symmetry breaking solution, 
becomes the probability distribution of these local field probabilities |M. This is a functional order 
parameter which is difficult to handle. Several interesting attempts at solving the one step RSB equations 
have been done in the past |3|,|2^,^, but they all restricted this fmictional order parameter to some 
particular subspace, within which a variational aproach was used. 

In this paper we present an improved solution of the Bethe lattice spin glass. This solution is nothing 
but the application to this problem of the cavity method, treated at a level which is equivalent to the one 
step RSB solution. It is valid for any connectivity and any temperature. In the next section we discuss the 
various tree-like lattices which are usually studied and precise our definition of the Bethe lattice problem. 
In se ct. PI we recall the basic steps of the simple Bethe-Peierls approach, and we discuss its instability in 
sect. |IV|. In sect. ^ where we discuss the formalism of the cavity approach at the one step RSB level. Sect. 
[Vl| describes the algorithm used to determine the d istri bution of local fields within this approach. The 



implementation of the algorithm is discussed in sect. VIl , where we derive explicit results for a lattice with 



six n eighb ours per point and compare the analytic prediction to those of numerical simulations. Finally, 



sect. VIII contains a brief discussion and mentions the perspectives. 



II. THE BETHE LATTICE 

We consider a system of N Ising spins, Ci ~ ±1, i G {1, ..., N}, interacting with random couplings, the 
energy being: 

E = - Y^ J^j(J,aJ . (1) 

The sum is over all links of a lattice. For each link < ij > the coupling Jij is an independent random 
variable chosen with the same probability distribution P{J). The various types of tree-like lattices which 
have been considered are: 

• A) The Cayley tree: starting from a central site i — 0, one builds a first shell of fc + 1 neighbours. 
Then each of the first shell spins is connected to k new neighbours in the second shell etc... until 
one reaches the L'th shell which is the boundary. There is no overlap among the new neighbours, 
so that the graph is a tree. 

• B) The random graph with fluctuating connectivity: for each pair of indices (ij), a link is present 
with probability c/N and absent with probability 1 — c/N. The number of links connected to a 
point is a random variable with a Poisson distribution, its mean being equal to c. 

• C) The random graph with fixed connectivity, equal to fc -f- 1. The space of allowed graphs are all 
graphs such that the number of links connected to each point is equal to fc -I- 1. The simplest choice, 
which we adopt here, is the case where every such graph has the same probability. 

On a Cayley tree a finite fraction of the total number of spins lie on the boundary. The Cayley tree is 
thus a strongly inhomogeneous system, the properties of which are often remote from those of a usual finite 
dimensional problem. For this reason people generally consider instead a Bethe lattice, which consists of a 
subset of the Cayley tree containing the first L' shells. Taking the limits L — > c», L' ^ oo with L/L' — > cxd 
allows to isolate the central part of the tree, away from the boundary. This procedure is OK when one 
considers a ferromagnetic problem. In the case of a spin glass this definition of the Bethe lattice is not free 
from ambiguities: one can not totally forget the boundary conditions which are imposed on the boundary 
of the Cayley tree, since they are fixing the degree of frustration p2]. For this reason we prefer to define 
the Bethe lattice as the random lattice with fixed connectivity (lattice C defined above). Clearly on such 
a graph the local structure is that of a tree with a fixed branching ratio. Small loops are rare, the typical 
size of a loop is of order log N. Therefore in the large N limit the random graph with fixed connectivity 
provides a well defined realisation of a Bethe lattice, i.e. a statistically homogeneous, locally tree-like 
structure. This is the lattice which we study in this paper (the case of fluctuating connectivities will be 
studied in a forthcoming work) . Numerical simulations of this system can be found in ]ll|,p2|,p8|pS|] . 



Historically, spin glasses on the Bethe lattice and on diluted lattices with a fixed finite connectivity 
(type C) were often discussed as a separate issue. The reason for these separate discussions of the same 
problem is the type of techniques which are used. Generally speaking the Bethe lattice papers rely on 
the Bethe Peierls method while the random lattice papers use the replica method. One exception is the 
use of the cavity method for the random lattice case p5| , p6[ . Hereafter we shall basically develop the 
iterative/cavity approach, but we shall also mention at each step its connexions to the replica approach. 

III. THE SIMPLE BETHE-PEIERLS 'SOLUTION' 

This section will give a brief review to the standard approach to the spin glass on the Bethe lattice, 
defined as lattice C in the above classification. This solution is wrong, because, as we shall see later, it 
does not consider the phenomenon of replica symmetry breaking, however it sets the stage for the correct 
solution that will be presented in the next section. 

A. The iterative approach 

As is well known, on tree-like structures the problem can be solved by iteration. Let us consider in 
general the merging of k branches of a tree onto one site ctq as in fig.0. The partition function can be 
computed exactly if one introduces, for each of the outside spins (Ji,i G {1, ...,fc}, the effective "cavity" 
field hi representing the action onto the spin ai of all the other spins, in the absence of the central spin ctq. 
In other words the magnetization of the i-th spin in absence of the central spin is given by mi — tanh(/3/ii). 
The variables hi (and consequently the variables mi) are uncorrelated in the limit N —* oo. It is crucial 
to consider the magnetization before the introduction of the spin ctq, because after its introduction all the 
spins that are coupled to the spin ctq become correlated. 
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FIG. 1. The merging of k branches off the tree (here fc = 3) onto the spin cro. The cavity field hi is the total 
field acting on spin ca in the absence of the central spin ao. 

Calling Ji the coupling between spins ao and Oi , the partition function of the spin ctq is expressed as 

/ k k \ 



y^ exp (3aQ ^ JiCJi + /3 ^ hiiJi 



(2) 



Cro,fTl,...fTfc 



Let us recall here the basic identity which allows to forward the effect of the fields hi onto spin ctq: and 
which is used repeatedly in this work. For an Ising spin ctq — il, one has: 



2_] exp (/3(To Jcr + Phcr) — c{J, h) exp(/3u( J, h)(7o) (3) 

cr=±l 

where we define the two functions u and c as: 

u(J,/i) = iatanh[tanh(/3J)tanh(/3/i)] ; c(J, /i) = 2 ^°'^^/^ ''°f/(^^^ (4) 

/? cosh[pu[J,h)) 

The magnetization on site is thus toq =< ctq >= tanh(/3/io), where 

fe 

/lo = ^u(-'^i,^j) ; (5) 

1=1 

From this equation one gets the basic recursion relation for the probability density Q{h) of local fields: 

k / k \ 

Q{h) = Ej ffl [dh.Qih,)] s(h-J2 u{J^, hi) j . (6) 

Here and throughout the paper, we denote by Ej the expectation value with respect to all the exchange 
couphng constants Jc Ej — ^Y\i [dhiP[Ji)\ . 

It will be useful for future use to introduce the probability distribution Riu) of the propagated field 
variable u{J,h): 

R{u) ^Ej I Q{h)dh 5{u - u{J, hj) (7) 

The field distribution is nothing but the convolution: Q{h) — J dui . . . duk R{ui) . . . R{uk)S{ui + . . . + 
Uk - h). 

In order to relate the true distribution of local fields, Qt{H), to the distribution Q{h) of local fields on 
one branch^ one needs to consider the merging of fc + 1 branches onto one site. The true local field Hq 
on a given site is simply given by a sum of contributions from each of its fc + 1 neighbours, 

fc+i 
Ho = Y.u{J„h,), (8) 

i=i 

where as before hj is the local field on j in the absence of the spin sq . This gives the distribution of true 
local fields Qt{H) as the convolution: 



, fe+i 

E 



Qt{H) = Jl[[du,R{^,)]5(H 



k+l 

Ui 
i=l 



(9) 



Let us now compute the internal energy with this method. We add a new link [Q with a coupling 
constant J^ between two spins ai and aj, where the local fields in the absence of the new link are 

respectively h^ and h): . Then the energy of this link is: 

Eij = -Jrj{(Tiaj) , (10) 

where the expectation value is computed using the Hamiltonian Hij(ai, <7j), which is given by 

H,j {a,,aj) = - [j^jo.aj + h^pa, + /if a^) . (11) 

A simple computation shows that 



'^We denote by upper case letters the true local fields and by small case letters the local fields on one branch. 



tanh(/3 Jy ) tanh(/3/ip^ ) tanh(/3/i['^ ) ' 



«i ~ "^^J . , .„^ . , .^.CiV . ,„,(i). ■ ' ^ 



Computing the total free energy of the system is sHghtly more involved. Using the fact that the Bethe- 
Peierls approximation is exact on the tree-like lattices one can write the free energy as the sum of site 
and bond contributions |14|-|l6||: 

^ = -^E^''+E^k' (13) 

where the contribution from the bond ij is 

- PF<%> = In E e^P i-PH.ji'^^.'^j)) (14) 

and that from the site i is: 

- /3i^/'^ = InE exp iPH.a,) , (15) 

where Hi is the total spin acting on spin ai. One can prove the validity of the expression ( |l3| ) by the 
following two steps: 1) it clearly gives the correct free energy at high temperature; 2) using the fact that 

X^iYi) ^i — kHi, where the sum is over all the neighbours j of site i, one finds that d{l3F)/dp gives back 
the correct expression for the internal energy obtained in (|l3) . We notice en passant that this free energy 
is nothing but the generalization to a finite coordination number of the TAP free energy (and reduces to 
the usual TAP free energy in the limit of infinite coordination number) [fsUlfl] . 

The Edwards- Anderson order parameter [Q, q = (1/iV) ^^ < at >^ can be written as the magnetiza- 
tion squared of a spin coupled to A: -I- 1 neighbours and is given by 

q= I dHQt{H) [tanh^(/3H)] . (16) 

We shall also compute the link overlap, g*-'' = {2/{N{k + l))X]<i7> "^ '^^'^i ^^' which is deduced from 
the Q{h) distribution as: 

gW =Ejf dhdh'Q{h)Q{h') panh(/3J)+tanh(/3fe)tanh(/3feO\ ^ 

^ -^ J ^V y^V '' l^i + tanh(/3J)tanh(/3/i)tanh(/3/(/)y ^ ' 



B. A variational formulation 

We have just seen in the previous section that, if one neglects the possibility of RSB, all the thermody- 
namic quantities of the Bethe lattice spin glass can be computed in terms of the probability distribution 
Q{h) of the effective field h. This probability distribution is obtained by solving the self-consistency 
equation (^. 

It is interesting for many reasons, some of which will become clear later, to reformulate this problem in 
a variational way. One can write a free energy i^[(3], which is a functional of the probability distribution 
Q{h), such that: 

1. The equation SF/6Q{h) = is equivalent to the self-consistency equation (|g) for Q{h). 

2. Calling Q* the solution of the previous equation, the equilibrium free energy m3t) is equal to F[(3*]. 
This free energy functional is given by: 

F[Q] - . ^ -• fc 



N 2 , 

i=l 

. k+1 



- / n [dh^dg^Q{h^)Qig^)] F^^) {hi . . .hk,gi . . . gt) 



k ffl [dh,Q{h,)] F'^'\hi . . . hk+i) (18) 



where 



PF^'\hi...hk+i)^Ejln 



■fe+i 



PF'-^\hi...hk,9i...gk)=EjEKln 



fJ^d{J,,h,) 



n 



J 0'Q,0'l,...CTfc_|_i 
1 



fc+1 



fe+l 



y^ exp Pap y^ J,fTi + /3 y^ ^i^i 

L i=l i=l 

E E 



cro,0"i,...crfc To,Ti,...rfc 



exp 



^■^ d{Ji,hi)d{Ki,gi 
/3 JoCToTo + /3cro ^ JiO-j + /? ^ h^(Ti + (3to ^ X^ri + P^ 9in 



(19) 



These two expressions are represented pictorially in fig. 0. In this formula the function d{J, h) is an 
arbitrary positive function, since its contributions to the two pieces F^^^ and F'^^ cancel. We shall mainly 
use it with d{J,h) = c(J,h), which allows an easy connection with the expression (13), but some other 
choice will also be useful in order to make contact with the result of the replica method, as we shall see 
below. 
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FIG. 2. A pictorial representation of the two contributions (h9[) to the free energy. The 'site' contribution on 
the top is obtained by merging fc + 1 lines onto one site (here fc = 3) , and the 'bond' contribution pictured on the 
bottom figure is obtained by adding one new link Jo, and two new spins ctq and tq, to the lattice, together with 
the other k branches arriving onto each of these spins. 

Let us now show that this free energy has the desired properties. In order to check that Q*{h) is a 
stationarity point of the free energy in the space of normalized probability distributions Q{h)[ such that 
J dhQ{h) — 1), we need to show that 6F/SQ(h) = constant when Q — Q*. This functional derivative is 
equal to 

h^§} " ^^^ + ^) I dh2... dhkQ{h2) . . . Q{hk) 

(J dgi... dgkQigi) . . . Q{gk) F^^\hi ...hk,gi... gu) - J dhk+iQ{hk+i)F^^^ {h, h2 

Using (||), one easily sees that, if Q{h) ~ Q*{h) satisfies the self consistency equation (g), one has, for any 
hi . . .hk-. 



l-k+l. 



(20) 



dgi... dgkQ*{gi) . . . Q*{gk)F^^^ {hi . . .hk,gi . . . gk) = A+ dhk+iQ*{hk+i)F^^^ {hi... hk+i) , (21) 
where A is a constant (independent of hi ... hk), given by 

A={-1/I3) fdgoQ*{go)E.,{ln[d{Jo,go)]+kln[c{Jo,go)/d{Jo,go)]) ■ (22) 



This shows that the functional derivative ( |20| ) is a constant. The repeated use of (||) allows to show 
similarly that the saddle point free energy F[Q*{h)] is indeed equal to the free energy (Oh (the factor 
{k + l)/2 in dig) is nothing but the number of links per site). 



As an extra check, one can see that the derivative of f3F[Q] we respect to /? gives the internal energy 
of the previous section (the derivative is quite simple if we absorb most of the /3's redefining the h and 
notice that the only explicit dependence on /3 comes from the term PJ). 

It is interesting to note that, using the basic recursion relation (|^), and the special choice d{J,h) — 
2cosh(/?/i), we can write the free energy under the simple form: 



F[Q] 



N 



- Y[[dh,Q{h^)] F^^'\hi...hk)-^- / l[[dh^Q{h,)] F^^\hi . . . hk+i) 



where 



■PF^^'\hi 



Ej In 



n 



2cosh(/3/i,) 



H exp 



CTo,a-i,...CTfc 



(3ao ^ JiCTi + /? ^ hia^ 



i=l 



i=l 



One must keep in mind that this expression for the free energy is correct only if Q satisfies eq. 
should not be used as a variational free energy (see also next section). 



(23) 

(24) 
and 



C. Equivalence with the replica formalism 



We shall not present here the details of the replica approach to this problem, for which we refer the 
reader to |l^,^,^,^,^ . Let us recall the main results of Q . In the replica approach one introduces a 
probability distribution p{a) , where the variables a are n Ising variables {n eventually goes to zero) . One 
can introduce a free energy functional Freplpicr)]. The equilibrium free energy is given by F{/3) — Frep[p*], 
where p* is the solution of the stationarity equation SFrep/Sp = 0. 

The expression of the free energy functional in replica space can be derived following exactly the same 
steps as in the previous section. The final result is p3| 



/3 rep 



[P] 



1 



■ In Tr„ 



p(cr)>(T)'=exp ^PJaaTa 



fcln(Tr, [p(a)'=+i]) (25) 



where Tr^. denotes the average over the 2" configurations of the variables a or t, and the correct result 
is obtained in the n — » limit. The same value for the free energy is obtained if we multiply the function 
/9 by a constant, so that the p does not need to be normalized, although it is more convenient to work 
with a normalized p. The advantage of the replica approach is that the system is homogenous and the 
distribution p{a) is the same in all the points. This advantage is partially compensated by the fact the 
the number of variables n is going to zero. 

As can be readily checked, p satisfies a very simple equation: 



P(^) 



EjTVr [p(T)fc exp (ELi PJaaTa)] 

Tr. [p{T)k] 



Using this equation the free energy (Eq) can be simplified to 



/3n 



-^ rep YP\ 

N 



^ln(Tr. [p[af]) - ^ In (Tr. Wf^']) 



(26) 



(27) 



where as before this new form of the free energy cannot be used in a variational formulation. The result 
( psj ) for the replicated free energy is correct in general, whether the replica symmetry is broken or not. 
The problem is to find the solution p*{cr). In the replica symmetric situation this task is easy: the p{a) 
is a function of only E = ^^ aa- We can thus write in general: 



p{a) — du R{u) exp(/3uS) . 
where the normalization condition of p{a) imposes: 

duR{u){2cosh{Pu)f = 1 



(28) 



(29) 



Using this expression for p, we obtain in the smah n hmit: 



^ fc+i 
In (Tr, [p{af+^]) ^ n \[ [du,R{u,)] In 

•^ 4=1 



■fc+l 

n 



'-J- 2cosh(/3ui) 



^ exp /3cro X! ' 



(30) 



and: 



In Tia.T 



In 



p{af pirf eyi^{ ^ PJcJaU) 



]^ [duidviR{ui)R{vi 



4=1 



n - 

^J- 4 cosh(/3ui) cosh(/3ui) 



y^ exp /3cro X! "* + ^'^" X! ""* + /3^oCTo 



To 



(31) 



Putting these expressions back into the rephca free energy (25), one gets exactly the functional F[Q] 
which we had written previously in (|l8|), provided we identify the n ^ limit of R{u) with the probability 
distribution (|^) of the variable m(J, h), and we use in ( p!^ ) a function d{J, h) — c{J, h)[2 cosh.{(3u{J, h))]. 

In other words we have seen three equivalent ways to solve the Bethe lattice spin glass in the replica 
symmetric approximation: 

• One can derive the recursion equations (pt) for the probability distribution Q*{h) of the local 'cavity' 
field h, and evaluate the free energy and the internal energy using this distribution. 

• Alternatively one can introduce the free energy functional (O), which depends on the probability 
distribution Q{h) and satisfies a variational principle: the distribution Q*{h) is obtained as the one 
which makes the functional stationnary. 

• One can obtain the same functional starting from the replica approach (P5|), making explicitely an 
assumption of replica symmetry, and doing some simple algebra. It is typical of the replica approach 
that a probability distribution is traded with a function of n variables, in the n ^ limit. 

D. Free energy shifts 

It may be instructive to compare this approach with the more usual cavity method and to check that 
we obtain the same results. In the cavity method, one computes the free energy by averaging the various 
free energy shifts obtained when adding a new site or a new bond to the lattice. 

The first quantity which we compute is the free energy shift AFner obtained my adding a new spin cfq 
connected to k branches, as we do in the iterative procedure. Using the same notations as in (ra), this free 
energy shift is: 



f3AFiter{Jl ■ ■ ■ Jk,hi . . .hk 



In 



2 cosh /3yju(Ji, hi) 



El- 



cosh(/3Ji) 



cosh(/3u(Ji,/ij)) 



(32) 



In order to compute the total free energy, we also need the free energy shift when adding the new spin 
CToj onto which merge k + 1 branches (see fig. 0). This free energy shift is equal to the same quantity with 
k changed into k + I: 



k+l 



/3AF(i) ( Ji . . . Jk+i ,hi... hk+i) = -/3F(i) {hi... hk+i) + PY. 1^[2 coshiPh,)] 



In 



fe+i 



2 cosh (iy^^u{Ji,hi) 



k+l 



cosh(/3Ji) 
cosh(/3u( Ji,/ij)) 



(33) 



The free energy shift when adding the two new spins cto, to (see fig.0) is equal to: 

k 

AF^^\ji ...Jk,Ki...Kk,hi... hk, 91... gk) = F^^\hi ...hk,gi... gk) ~ Y. ln[4cosh(/3/i,) cosh(/3g. 



(34) 



and is given by: 

A; 

- /? A F(2) ( Jo, Ji . . . Jfc, Xi . . . Xfc , /ii . . . /ifc , 51 . . . fffc) = ^ In 



In 



4=1 

k 



cosh{(3Ji) cosh{f3Ki) 

cosh{Pu( Ji, hi)) cosh{l3u{Ki,gi)) 



(k k N 

PJoo-QTo + (3aQ ^ u( Ji, hi) + /3ro ^ u{Ki,gi) 
1=1 j=i / 



(35) 



In the process of adding new sites or new bonds randomly, one can thus compute the total free energy 
as the average over the distribution of fields and couplings of [(fc + l)/2] Ai^^^^ — kAF^^\ It is a simple 
exercise to check that this indeed gives back the free energy (118). 

IV. THE RSB INSTABILITY 

The recursion relation of the local fields (^ has been the subject of a lot of studies in the past twenty 
years p^- p^ , ^ , pl[ . The distribution Q{h) is a simple S function at the origin, indicating a paramagnetic 
phase, at high temperatures (3 < (3c, where the critical inverse temperature /3c is the solution of |17t : 

Ejtanh^{PcJ)^l/k . (36) 

In the low temperature phase the specific heat becomes negative at low enough temperatures at least for 
some distributions of couplings |1^, and the solution for Q[h) becomes identical to the replica symmetric 
field distribution of the SK model in the large k limit |1J,0,0| , which is known to be wrong. Another 
indication that the above procedure gives a wrong result for the Bethe lattice (while it might be correct 



for the Cayley tree with some sets of boundary conditions |18 2lf|) is the fact that it fails to identify a 
transition in a magnetic field H. This transition exists though, on a line in the H — T plane similar to 
the A-T line, and can be identified by considering the onset of correlations between two replicas of the 
system 0. 

One can investigate the instability of the previous solution using the replica method. Writing the 
recursion relations for the replicated system, Mottishaw has shown that the replica symmetric solution, 
which coincides with the simple Bethe-Peierls iteration described above, is unstable at (3 > f3c (or, in a 
field, beyond the A-T line) |l3]. Unfortunately, getting the replica symmetry broken solution in the low 
temperature phase is difficult. In general the problem involves an infinity of order parameters which are 
multi-spin overlaps |24J3,n3,E0,p3[ . As we saw, the replica symmetric solution already involves an order 
parameter which is a whole function (the distribution of local fields); going to a 'one step RSB' solution 
P|, the replica order parameter becomes now a functional, the probability distribution over the space 
of local field distributions [g^ (the reason will be discussed in details in the next section). While one 
can write formally some integral equations satisfied by this order parameter, solving them is in general a 
formidable task. The solution is known only in the neighborhood of the critical temperature |22l, or in 
the limit of large connectivities p3| , where the overlaps involving three spins or more are small and can 
be treated perturbatively, allowing for some expansion around the SK solution. In the general case, the 
only tractable method so far has been an approximation which parametrizes the functional by a small 
enough number of parameters and optimizes the one step RSB free energy inside this subspace p3,p3] . 
We shall develop in the next two sections a solution to this problem. 



V. THE FORMULATION OF THE 'ONE STEP RSB' SOLUTION 

A. The iterative approach 

In this section we shall explain the physical nature of the RSB instability and work out the equivalent 
of the 'one step RSB solution' using the cavity method [Hj3], i.e. the same type of iterative approach 



which was used in sect. [II A 



The reason for the failure of the simple iterative solution is that it neglects the possibility of the existence 
of several pure states H . We shall proceed by first assuming some properties of the states on one branch, 
and then imposing the self-consistency of these hypotheses when one joins k branches to a new site. Let 



us assume that there exist many pure states, labelled by an index a going from 1 to oo, with the following 
properties: Looking at one branch of the tree as in fig.l, the total local field seen by the site i = at 
the extremity of this branch depends on the state a and is denoted by h^ ; the free energies of the states 
on one branch, i^", are independent identically distributed (iid) random variables, with an exponential 
density behaving as 

p{F) = exp(/?a;(F - ^^)) , (37) 

where F^ is a reference free energy. The differences between the free energies remain finite when the 
volume goes to infinity, which means that the various states have non-zero statistical weights: 

E,exp(-/3F'r) ' ^'"^ 

The fact that the W can be normalized in this way is possible only if a; < 1 . 

Let us consider as before a point i on one branch of the Bethe lattice, i.e. a point connected to k 
other points. In each phase a of the system the magnetization to" will be different and therefore the 
effective fields hf^ depend on a. The description of the properties of this point will include the list of 
fields /if and the free energies of the branch F", for all states a. Here we shall assume a relatively simple 
situation namely that the free energies and the magnetic fields are not correlated, and the distribution of 
free energies is the one described above, leading to the density (p7|). It is convenient (to avoid possible 
difhculties in dealing with measures in infinite dimensional spaces) to order the states in an increasing 
order of free energy, and to consider only the set of the first M. states with lowest free energies (in the 
end M. will be sent to infinity). Let us introduce the set of the local fields in all the M. states, h = {/if}. 
When changing the sample (or equivalently changing the site i), these fields fluctuate and our task is to 
compute the corresponding probability distribution (3(h), which is a function of the M. effective magnetic 
fields which is left invariant by the permutations of these fields. This task is simplified if we assume that 
theA^ fields /if on one point can be characterized as independent random variables. We thus assume that 
there exists a probability function Qi{\\) which can be written in a factorized form: 

M 
Qi(h) = [] 0,(/i„) . (39) 

a = \ 

The total probability function Q(h) is thus given by 



N 



Q(h) = iV-i^ 



M 



n Q«(^") 



(40) 



In other words on any given point the fields are independent variables, which become correlated on the 
global level after we average over the samples. More generally one can assume that the total probability 
function (3(h) is of the form: 



„ M 

Q(h)= /dAm(A) JJ(7(/i"|A) 



(41) 



where A is an appropriate set, to(A) is a probability distribution, and q{h\\) is a probability distribution 
on /i, conditionned to a given value of A. A possible representation of the distribution ( [4l| ) is given by 
(pO|), where each point i is characterized by a parameter A (extracted with the measure r7i(A)). 

We shall now check that this hypothesis is self consistent, i.e. that it is reproduced when one iterates 
the construction of the tree by merging k lines to a new spin a^ [p5|. For each state a, the local field on 
this site /iq is expressed in terms of those on the branches, hf by^5|), giving: 

k 

/i« = ^zi(J„/jf). (42) 

The free energy shift AF" for the state a during this process is given by the function ^Fuer defined in 
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AF'^ = AF,ter(Jl ... Jfe, /!?,... , K) . (43) 

We must be careful at this stage because the free energy shifts and the local fields on the new spin a^ 
are correlated. More precisely, for a given state a, h^ and AF" are two correlated variables, but they are 
not correlated with the local fields or free-energy shifts in the other states. Because of our ordering process 
of the free energies, we need to compute and order the new free energies G" = F" + AF". G" and C 
in two different states are obviously independent random variables. Furthermore, a standard argument of 
the cavity method 0,0], relying on the exponential distribution of the free energies, allows to show that 
the new free energy G" is in fact uncorrelated with the local field Hq. To show this, let us introduce the 
joint distribution Po{ho, AF) of Hq and AF". The joint distribution Ro{hQ, G) of the local field and the 
new free energy is given by: 

Roiho, G) oc / dFd{AF) exp(/3a;(F - F^j)Po{ho, AF)S{G - F - AF) oc exp{f3x{G - F^))go(/io) , (44) 

where 

Qoiho) = C j d{AF) Poiho,AF)eM-0xAF) (45) 

the constant G being fixed in such a way that Qo{ho) is a normalized probability distribution. 

In our ordering process of the new free energies Ga we pick up the A4 lowest ones, sending in the 
end A4 to infinity. This ordering process thus gives rise to a probability distribution for the Hq in which 
the fields for different a are not correlated and have the distribution Qo{hQ). The reader should notice 
that this distribution is in general different from the naive result JdAFg Po{h, AFq). In this way we 
have constructed, for one given new spin ctq with a fixed environment of coupling constants, the new 
distribution of all local fields: Yla Qoihg)- By averaging over the coupling constants, one thus generates 
the probability distribution Qo{h.) which is a functional of the probability distribution Q(h) of the other 
k sites. Imposing that 

Qo(h) = Q(h) (46) 

gives a self-consistency equation for the probability distribution (5(h). We shall see in sect. W^ how one 
can actually find a solution to this self-consistency equation with a good accuracy. 

Let us suppose for the time being that we know the self-consistent distribution Q(h) and let us evaluate 
the free energy and the internal energy. The computation is quite similar to the one in the replica 
symmetric case. There are two contributions to the free energy: the site contribution and the bond 
contribution. 

The local site contribution to the free energy is evaluated as a weighted average of the free energy shift 
when adding one new spin: 

F(i) = -^Ej{\n f £ W^ exp[-/3AF(i)(Ji . . . J^+i, /i? . . . h^+,)]\) , (47) 



where AF*^^^ is the site free-energy-shift function defined in (pSh, and the bracket denotes an average over 
the distribution of weights W" derived from ( [3^j37| ), and over the fields (with distribution ni=i QO^i))- 
The local bond contribution to the free energy is evaluated as a weighted average of the free energy 
shift when adding a new bond and the correspoonding two spins: 



^EjEK{\n I f2 W' exp[-/3AF(2)( Jo, J, . . . Jk,K, . . . Kk,h'^ . . . K, g^ . . . g'^)]\] 



where AF^^^ is the bond free-energy-shift function defined in (p5|), and the bracket denotes an average 
over the weights, and over the fields (with distributi on F fi-i [Q(hj)Q(gi)]). 
The total free energy is given, according to ( |l3| , |l5| , |l4|) , lay: 



F = :^±iF(2) - A;F(i) . (49) 
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Similarly to what happened in the RS case (^3|), one can also write here a simplified form of the free 
energy, valid only on the saddle point: 

F = :^±1f(i') - :^^f(i) (50) 



where 



1 f ^ 



■K)]])- (51) 



For reasons which are beyond our control this second form of the free energy has empirically smaller errors 
(about by a factor 3) and less systematic errors (at finite Ai) than the first one. 

The computation of the internal energy is done by considering what happens when we couple two sites 
which wherepreviously connected each to k branches of the tree. We obtain, with the same notations as 
in formula (Qq) and in fig. 0: 



M 



tanh(/3 Jo) + tauh{f3h^) tanh(/35^) 



t^i ^ "l + tanh(/3Jo)tanh(/3/i^)tanh(/35^)^ ^ ' 

where /ig = X]i=i "^{Ji-, ^f )? 5o ~ X]i=i "(^ij 5f ) ^^^ we have introduced the shorthand notation: 

^. ^ P^°exp[-/3Af(2)(Jo,Jl...Jfe,gl...gfc,/^^..feg,g^..gg)] 
' E^iW^''exph/3Ai^(2)(j^^j^...j^^^^...^^^;,7,.j,7^^7,,,^7)] 

The various overlaps can be obtained in the same way. There are now two site overlaps, the self-overlap 
gi and the inter-state-overlap go? which are given by: 

M fe+l 

a— 1 z— 1 

go = Ej^Y. VFfM/f tanh[/3^M(J„/in]tanh[/3^u(J„/if)]) (54) 

where we have kept the same notations as in (|47| ) but we have introduced W^ which is a notation for: 

^„ ^ exp[-/3f " - /3Af (1) ( Ji . . . Jfe+i, /.? . . . feg^,)] 
' "E^exp[-/3i^7_/3AF(l)(Jl...J,+l,/^^../^^+l)] • ^ ^ 

Similarly, we have two link overlaps g)^ and gp which are given by: 



ID ^ p / V V7" / tanh(/3Jo)+tanh(/?/ig)tanh(/3gg) \ \ 

^1 '^^^ 2 ll + tanh(/3Jo)tanh(/3/i")tanh(/3g?)/ ' ^ ' 



tanh(/3Jo) tanh(/3/i^) tanh(/35^) 
and 



(0 _ z^./Y- w«w/3 /^ ta nh(/?Jo) + tanh(/3fag) tanh(/?gg) 

tanh(/3Jo) tanh(/3/ig^) tanh(/3g^) 






tanh(/3Jo) -h tanh(/3/i^) tanh(/3g^) \ 



1 + tanh(/3Jo) tanh(/3/i[^) tanh(/3g^) 



where we keep the same notations as in (B2Lp3[). 

At this stage we have written the hole formalism of the cavity method at the level of one step RSB. 
The self-consistency equation (p6| ) fixes the distribution Q(h), from which one can deduce the free energy 
and internal energy through ([49|) and ( p2[ ) . One can actually find a self-consistent solution for any value 
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of the parameter x £ [0, 1], and the free energy depends on x through the distribution of free energies. It 
is well known that, in order to describe the thermal equilibrium, one must fix x by maximising the free 
energy with respect to x |||,^,Q . (Actually the whole dependance on x carries some information p6|~p8| , 
particularly interesting for optimization problems, which we shall not try to study here since this paper is 
restricted to the study of equilibrium thermodynamics). As we shall see, the variation of free energy with 
respect to x is small and one needs a better computation of the derivative than just a naive difference. 

We have improved the precision on the computation of the free energy and its x derivative by using the 
theorem of appendix 1 , which allows to compute explicitely the derivative of the free energy with respect 
to x. From the structure of eq. (M9) one finds that the total derivative d — dF/dx takes the form : 



d(x) = --i^-fcd(i) + ^d(2) 



(58) 



where 



d(i) = -^Ej{\n If^W' exp[-/3AF(i)(Ji . . . J,+i, /^^ . . ht+,)] HAF^'\j, . . . J^+uK . . . K^,)]\) 



(59) 



using the notations of (|47|), and 



rf^'^ ^.-'^EjEK{\n f £ M/"exp[-/5Ai^(2)(Jo, Ji ...Jt.K,. ..Kk,h^ ■■■K,9?- . . g^)]HAF('^]\) , 



(60) 



using the notations of (Eq). 



B. A variational formulation 



As in the case where no replica symmetry breaking is present we can write a free energy functional of the 
field distribution Q(h) such that the self-consistency equations for Q are equivalent to the stationnarity 
condition of this functional. 

This free energy functional is a simple generalisation of the replica symmetric one, given by: 



F[Q] 
N 



k + l 



J l[[dh,ds^Q{h,)Q{g^)] F(2)(hi . . .hfe,gi . . .g,) 



I, k+l 

kj Y[[dh,Qih,)] FW(hi...hfc+i) 



(61) 



where 



/?i^« (hi . . . hfc+i ) = £;j (In K] W^" 



■fc+1 

n 



^\ 2cosh{phf) 



/3F(2)(hi . . . hfc, gi. . . gfc) = EjEK{\n Y. W 



n 



Y. ^^p 



ao,cri,...(Tk^i 



fc+1 



fc+1 



/3(7o^J.a,+/3^/ifa, 



i=l 



i=l 



), 



(62) 



E E 



CTO:""! ■•••<'■ Ic TO.Tl.-.-Tfc 



exp 



L 4cosh(/3/if)cosh(/55f) 
f3JoaoTo + (3ao ^ Ji(J^ + (3^, K^^i + /?^o E ^*^* + /^ E df-^^ 



I). (63) 



The weights Wa are given by (pq), and the brackets stand for an average over the distribution of free 
energies (p7|). 
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The proof of the equivalence of the stationarity equation of this functional with the self-consistency 
condition of the iterative procedure can be done exactly as in the replica symmetric case. The advantage 
of this variational formulation is that we can compute directly the various derivatives of the free energy 
(e.g. with respect to x and with respect to [3) by taking into account only the explicit dependence. 



C. Equivalence with the replica formalism 



We can compare what we have obtained in the previous sections with the results from the replica 
formalism. In the one step RSB formalism the n replicas are divided into n/x groups (labeled by C) of x 
replicas each, and the function p{a) depends on the n/x 'block' variables 



^c = / , gg 1 



(64) 



a^C 



each sum containing x terms. 

In this case we face the problem that yo(cr) may depend on the n/x variables Sc in a rather complex 
way. Similarly to what we have done in the iterative approach, we shall not try to describe the most 
general dependence, but we restrict to the class of probability distributions p{a) which can be written as: 



/ n/x I n/x 

d\p{\) / W [dMc0(uc|A)]exp /3^ucS, 
•^ C=l \ C=l 



(65) 



with a positive probability distribution /i(A) ||32| and a positive function ip{u\X). The normalisation 
condition on p{a) is implemented by imposing that: 



VA 



so that the function 



du(j){u\X){2cosh{Pu)f = 1 , 



$(u|A) = (/)(M|A)(2cosh(/?M))" 



(66) 



(67) 



is a probability distribution on the u variable, for any value of A. 

It is easy to check that this form for the function is consistent with the stationarity equations for the 
free energy (Eq) by proving that, if p(o') has the form (pq), so does 



EjTiw p{Tf exp[^ fiJaaTa] 



(68) 



Using the Ansatz (pq) for p, one can write the 'site' term in the replica free energy (Ea) as: 



Tr. {p{at^') 



■fc+i 

W^ dXipiXi 

.1=1 



A(Ai,...,Afc+i)"/^ 



where: 



A(Ai,...,Afe+i) = / [|[dui$(u,|Ai)] 



2cosh(/3E-=i'^») 
nj2cosh/3u,] 



Using the small theorem proven in appendix I, the previous expression can be written as 



In^ (A, . . . , A.,0 - -/ n n M<<J>(<|AO] ( [in (^E ^"^f^|g^^ 



(69) 



(70) 



(71) 



where the bracket means an average over the distribution of the weights Wa have the usual distribution 
of the weights over states, parametrized by the parameter x (see appendix I). One gets in the end: 
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„fc+l 
InTr, {p{af+^) = " / H \-d\,p{X 



I n n [^<*("fiA.)i (.» (e »'"^ifS|?r) > '"' 

i— 1 a \ a i xi L I I A y 

The previous quantity is exactly the expression of the 'site contribution' to the variational free energy 
found in (|6^), provided we identify the probability distributions defined in the iterative approach (^l|) and 
in the replica approach (|65|): 

/i(A)=m(A) ; ^{u\X) ^ Ej f dh q{h\X)6{u - u{J,h)) . (73) 

A similar computation shows that the 'bond' term in the free energy, calculated either through the iterative 
procedure or through the replica method also coincide. 

We have thus derived the variational equations for the probability distribution of local fields in the 
one step RSB case using two different methods, the cavity iterative approach on the one hand, and the 
algebraic replica formalism on the other hand. 

VI. SOLVING THE ONE STEP RSB EQUATIONS 

Our method consists in following the population of local fields hf when one iterates the merging process 
of k branches onto one site. In some sense it thus amounts to solving the complicated equation for the 
functional order parameter by a method of population dynamics. In other words we parametrize the 
probability distribution by presenting a large number of instances of variables distributed according to 
this probability distribution. A similar method has ben first used in the context of mean field equations 
in ||]. 

To explain it in more details, let us first state how the procedure works in the case of the 'replica 



symmetric' approximation of sect. [II. There, one just chooses a population of TV local fields hi. At each 
iteration, one picks up k such fields at random among the TV, and computes the new field ho according 
to (ra). Then one field is removed at random from the population and is substituted by hg. 

In this way one defines a Markov chain on the space of the TV magnetic fields. This chain has a 
stationary distribution which is reached after some transient time. In the tV ^ oo limit, the stationary 
distribution satisfies the self-consistency equation(g). It is possible to argue that the corrections to this 
limit are proportional to N'~^ and could also be computed analytically. Our procedure consists in fixing 
the value of TV, iterating the merging transformation many times in such a way as to obtain the average 
over the asymptotic distribution at fixed TV with a high precision, and finally extrapolating the results to 
Af ^ oo. 

If we consider the case where there exist many states, we have the same problem as before, with the 
only difference that at each point we have a probability distribution Qi{h). We must therefore consider 
a population dynamics in which the elements of the population are probability distributions. We use the 
population method to represent the probability distribution in each point i by a populations of fields. In 
this way we have a population of TV populations of M fields (a total of AfM fields /if, i G {1, ...,TV}, 
a e {1,...,TV4}), where both TV and TW have to go to infinity. 

The basic step of the algorithm is the merging of k lines. One chooses k sites ii,...,ik in {1, ...,TV}, 
and one generates, for each of the Ai states, a new local field h^ obtained by merging k branches, using 
(||), as well as the corresponding free-energy shift AF" calculated using (|3|). But the field h^ is not the 
one which will enter the population of fields. The reason is that one needs to reweigh the various states 
by a factor which depends on their free energy shifts: as seen in (Ea), the field distribution, at a fixed 
new free energy, is modified by a factor exp(—f3xAF). From the knowledge of the hg one can infer an 
approximate form of the distribution Po^ho) from which they are extracted. For instance a simple form 
for Pq is a smoothly interpolated version of the identity 

/ Poiho)dho = (l/M) Y, ^(ho - K) (74) 
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where 9{x) is Heaviside's function (in practice we smooth this staircase function by a linear interpolation 
procedure). According to ([45|), the real field distribution Qo{h) is well approximated by a smoothly 
interpolated version of the identity 



/ Qa{h)dh = [l/M) ^ exp(-/3xAF")6i(/i - h^ 

J — OO 



(75) 



We can use two different methods in order to achieve the reweighing, which will lead to two different 
algorithms. 

• Method A: The idea is to generate, from the set of A4 fields h", on each of the sites ii, ..., ik which 
is used in the merging, a larger population, of rAi local fields hf , a' = 1, ...,rA4, taken from the 



same distribution. This will be realised by having 

M 



rM 



{l/M) J2 (^ih - K) ^ i^/irM)) J2 O^h - hi) , 



(76) 



at the level of linearly interpolated functions. Simultaneously one generates rM. independent random 
free energies F" , a' — f , ..., rM., with the exponential density (p7\). For each of the rM states, one 
then compute the local field hg and the free energy shift AF" . The correct reweighing is obtained 
by the selection of low lying states: one computes the new rM free energies F" + AF" , orders 
them, and keeps only the M states with the lowest new free energies. Their local fields, called h", 
with a € {1,---,A^}, have the correctly reweighed distribution provided that r is large enough so 
that there is a zero probability for the states a' with the highest free energy F" to enter the list of 
the M selected states after reweighing. 

• Method B: The idea is to generate directly the fields with the reweighed distribution ([Tq). Knowing 
the M local fields h^ and their free-energy shifts AF", we generate M new local fields h" in such 
a way that the following identity holds at the level of linear interpolation: 



J2a exp(-/3a;AF° 






exp{-pxAF°')9{h -K)-JlYl ^(^ 



M 



/i") . 



(77) 



Having generated the new fields h" which are typical of the properly reweighed distribution, we then 
substitute in the population the set of M local fields hf, a = 1,..,M by the set of new fields h". The 
site index i on which this substitution is performed is chosen sequentially. 

While the merging of k lines is enough to build up the Markov chain which generates the population of 
local fields, one also needs to consider some different merging processes in order to compute the various 
observables, free energy, energy, local overlap and link overlap. 

By merging fc + 1 lines instead of k, one generates with exactly the same procedure as above the three 
sets of M local fields /ig , free-energy shifts AF" and new fields ft," (we call new fields the fields which are 
typical of the reweighed distribution, obtained either through procedure A or B). Using the little theorem 
of Appendix I, the site contribution (Wm to the free energy is computed as 



F(i) = -(f/(/3a;))ln 



and the site overlaps receive a contribution 



(f/Al)^exp(-/3a;AF°) 



qi^Yl tanh2(/3/i") ; qo = Y^ tanh(/3/i") tanh(/?/i^) 



The X— derivative of the free energy (M) receives a site contribution equal to: 



S' 



-(l/(/3.x))ln 



(l/M) Y^ exp(-/3a;AF") ln[AF° 



(78) 



(79) 



(80) 
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These contributions are then averaged over many iterations. 

One can also add a new Unk to the system. This is done by merging k hnes on the left side of the 
new link, generating the local fields /iq , and similarly merging k lines on the right side of the new link, 
generating the local fields gQ h^. The corresponding free energy shift AF" is now computed using (p5|). 
The bond contribution (Esh to the free energy is computed by the same expression as (uq), but using this 
free energy shift AF"' corresponding to a bond-addition. 

The X— derivative of the free energy (pQ) receives a bond contribution (BG) equal to: 



d(2) = -(l/(/3a;))ln 



(l/M) V exp{~pxAF") ln[AF° 



a 



(81) 



In order to compute the internal energy and the link overlaps, it is useful to introduce the local fields Vq 
defined by 

„ ^ 1 , tanh(/3Jo) + tanh(/3feg) tanh(/3gg) 

^0 ^atan i ^ ^ ^^^^pJ^•^ tanh(/?/i«) tanh(/3<?^) ' ^ ' 

From this population one generates a population of new fields w" which takes into account the appropriate 
weights of the fields with one of the two procedures A or B, similarly to what was done in (JTq) or in ( [77[ ) 
for going from the fields h^ to h". From (^ 56|57| ), the contributions to the internal energy and link 



overlaps are given by: 

U ^-JaJ2 tanh(/3z;") , (83) 

a 

q[''> = ^tanh2(/3f") , g^'' = Y^ tanh(/3w")tanh(/3w'^) . (84) 

In order to compute the free energy with the simplified formula ( pO|) , one also needs the contribution F^^ 



to the free energy, which is obtained by the same expression as (78), but using the free energy shift AF"' 
obtained by merging k lines. 

Let us therefore summarize the main lines of our population dynamics algorithms for solving the Bethe 
lattice spin-glass problem at the level of one-step RSB. We have used two algorithms, A and B, which 
differ in the reweighing procedure used, but have otherwise the same skeleton: 

1. Start from the population of A/" x A^ local fields hf. 

2. Merge fc -I- 1 lines and compute the site observables: 

a) Choose at random fc -I- 1 sites «i, ...,ik+i in {1, ...,7V}. 

b) For each of these fc + 1 sites, say on site j G {ii, ...,ik+i}, one has a population of A4 local 
fields /i", a — l,...,A/(. For each of the A4 states, compute the new local field Hq obtained by 
merging fc -I- 1 branches, using (g). Compute the corresponding free energy shift AF" using (|33|). 

c) Knowing the sets of fields H^ and free energy shift AF" , generate a new set of fields H" 
according to ( [Tq ) in algorithm A (resp.(|77|) if one uses algorithm B). 

d) Compute the site contribution to the free energy using (|7q), its x— derivative using (jsy) and 
the contribution to the site overlaps (y%. 

3. Merge 2fc lines onto a new bond and compute the bond overlaps: 

a) Choose at random 2fc sites ii, ..., z^, ji, ..., j^ in {1, ..., A/"}. 

b) From the sites ii, ...,?fc, compute the A4 local fields h^ obtained by merging the fc branches, 
using M). From the sites ji,...,jk, compute the Ai local fields g^ obtained by merging the fc 
branches, using (g). Deduce the M local fields Vq using (^2|). 

c) Compute the free energy shifts AF°' using (|35|). 

d) Knowing the sets of fields Vq and free energy shifts AF", generate a new set of fields w" 
according to ( [76| ) in algorithm A (resp.(|77|) if one uses algorithm B). 

e) Compute the bond contribution to the free energy using (|78|) , its x— derivative using (rTJ) 
and the contribution to the internal energy ( p^ ) and link overlaps 
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4. Merge k lines and update the population of fields: 

a) Choose at random k sites zi, ...,ifc in {1, ...,Af}. 

b) For each of these k sites, say on site j £ {ii, ..., ik}, one has a population of A4 local fields 
h", a — 1, ...,A4. For each of the A4 states, compute the new local field Hq obtained by merging k 
branches, using (h). Compute the corresponding free energy shift AF" using (p2|). 

c) Knowing the sets of fields h^ and free energy shift AF", generate a new set of fields h"' 
according to (Um in algorithm A (resp.(|77|) if one uses algorithm B). 

d) Compute the contribution F'^^ ^ to the simple form (KG) of the free energy using (uSh. 

e) Pick up the site « £ {1, ..., A/"} sequentially, and substitute the local fields hj, ...., h^ by the 
new local fields ft,". 

5. Start again the iteration from 2). 

Obviously one needs not really perform the above three merging procedures sequentially. In our actual 
algorithm, we select 2fc + 2 random points, merge two groups of fc + 1 to compute site ovbservable, two 
groups of k to compute bond observables and to update two new sets of M. fields. 

A word about the difference between the two algorithms. In algorithm A, the value of r must be chosen 
large enough so that the probability of a state with the highest old free energy F" to enter the set of 
selected AA states be negligible. Both M and r must go to infinity and in the numerical computations 
we have taken r ~ M.. Algorithm B is faster than algorithm (A) by a factor that is asymptotically 
proportional to r for large r. In algorithm B there is no need of introducing rM. fields at an intermediate 
stage: it corresponds to the discretization of eq. ^ and we care take of the effects of the reshuffling by 
explicitly reweighing the fields. Unfortunately, as we shall see in the next section, the finite M corrections 
are empirically larger in algorithm B that in algorithm A. In our case algorithm B turned out to be 
faster by a factor about 10, but we mentioned both algorithms because algorithm A is somewhat simpler 
conceptually (and closer to the original discussion of the cavity method), and also because in different 
situations (e.g. depending on the value of x) the relative advantages may be reversed. 

VII. NUMERICAL RESULTS 

Here we present the numerical results for one case in order to study the dependence of the algorithm 
on the various parameters involved in the numerical computation. 

We consider the Ising spin glasses with binary couplings (J — ±1) on a random lattice of fixed coor- 
dination 6 (k=5) . High precisions measurements p9[ of the internal energy are available at temperatures 
greater or equal to .8 for different values of the number of spins N (up to A^ = 4096). The data of the 
energy at T = .8 versus size can be very well fitted by a power law correction: 

U{N) = U + A7V~" (85) 

with a quite reasonable value of w ~ .767 ± .008 and A ^ 2.59 ± .02. The value of the internal energy at 
infinite size is estimated to be 

[/ = -1.799±.001 . (86) 

We have done a replica symmetric computation for different values of the population size TV. For large J\f 
there are corrections proportional to 1/A^ (as expected) and for TV > 10^ the 1/TV corrections are negligible 
within our accuracy. With / — 100, S — 1000, TV — 4000, we obtain the following replica symmetric results 
for the free energy, internal energy, entropy, site and link Edwards- Anderson parameters: 

F = -1.863 ±0.002 , [/ = -1.8160 ±0.001 , 5 = .058 ± 0.004 

q = 0.6863 ± 0.0002 , qu^k = 0.6385 ± 0.0003 . (87) 

Notice that the value of the internal energy totally disagrees with the one found in the simulations (pq) . 
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3. The a;-derivative of the free energy, d{x), evaluated at a: = .21, is plotted as function of A^ ^ for the 
algorithm A (upper curve) and the algorithm B (lower curve). 



We have done a computation at the one step RSB level using the two algorithms described in the 
previous section, always at temperature T = .8. 

The crucial point is to find the value of the parameter x such that the derivative of the free energy 
with respect to x vanishes. In fig.|3| we show our results for the derivative d{x) at x = .21, plotted 
versus different values of A4. We have used both algorithms A and B. With algorithm B we have used 
A4 = 2',/ = 3 ... 12 and we plot the results obtained for / < 10 for a clearer figure. With algorithm 
A we have used r = Ai < 400. Unfortunately the finite A4 corrections are empirically much larger 
in the a-priori-faster algorithm B. Moreover, in this case although the finite A4 corrections seem to be 
asymptotically proportional to 1/A^, high order corrections cannot be neglected unless A4 is very large. 
In the end both algorithms give compatible asymptotic results at large Ai as seen on fig.pi with similar 
computer efforts (for this temperature and the values of x which are relevant). 

In fig.^ we plot the extrapolated result at A^ = cx) for the free energy derivative d{x), obtained using 
algorithm B. The data has been extrapolated with a second order polynomial of A^~^ in the interval 
M = [256-4096]. 
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FIG. 4. The a;-derivative of the free energy , d{x), computed using the algorithm (B) and extrapolated at large 
values of A4, plotted versus x. 

Replica symmetry breaking is clearly present. The value of x where the free energy is maximum, which 
can be obtained by estimating the zero of the function d{x), is given by x* = .24 ± .02. We use a similar 
procedure for all other observables: we extrapolate the x dependent results at A^ = oo and evaluate the 
errors due to the imprecise location of x* (which is by far the largest source of error). This gives the 
following values of the free energy, internal energy, entropy, site and link overlaps: 
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-1.858 ±.002 , U 



-1.799 ±0.001 , S' = .074 ± 0.004 
qi = 0.779 ± .006 , qo = .30 ± .02 



{1} 



.706 ± .007 , ql}^ = .408 ± .01 



(89) 



The energy is in very good agreement with the results from simulations. In order to compare the 
values of the overlaps, one can study the quantity (g^) = / dqP{q)q^. In our RSB theory we find (g^) = 
(1 — x*)qf + x*qQ — .485 ± .015 which agrees well with the numerical value (q^) — .49 ± .02. In order 
to perform a finer comparison it is useful to consider a quantity which is sensitive to replica symmetry 
breaking. A natural choice is 



R = 



j dqP{q)q^ - (j dqP{q)q' 



(90) 



We find in our RSB theoryi? = .046 ± .002 which is again in good agreement with the result of the 
simulations extrapolated at infinite volume: R — .051 ± .002. The agreement is remarkable if we recall 
that in the replica symmetric case R — Q. The possible small difference between our value and the 
simulation data is likely due to the the approximation of one step replica symmetry breaking (it is quite 
likely that replica symmetry should be broken in a continous way, as happens in the limit of infinite 
coordination number). One should notice that the order parameter q^ is non zero, which explains why 
some previous attempts at solving the one step RSB problem within a restricted subspace with go = did 
not improve much on the RS solution [|32| (the necessity of having a non vanishing go was already noticed 
in (§]). 

Finally let us point out the crucial effect of the reweighing of the states which modifies the flocal fields 
as in (^5|). In fig. (||) we plot the probability distribution of the field i?o before and after the reweighing, 
at a; = 1. 
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FIG. 5. The probability distribution of the field Ho at a; = 1 before (full curve) and after the reweighing (dashed 
curve) . 



VIII. CONCLUSION 

We have presented a general solution of finite connectivity spin glass problems at the level of one step 
RSB. This solution uses the cavity method, together with a kind of population dynamic algorithm to 
solve the complicated functional equations. As exemplified by a detailed numerical study of the sping 
glass on a random lattice with fixed connectivity, this method aloows to obtain good agreement with 
numerical simulations. It should be rather easily extendable to many other disordered systems with a 
finite connectivity, including the fluctuating valence spin glass and various optimization problems such as 
the K-satisfiability problem. 

The possibility to go to higher order of RSB should be explored. In principle the method we have 
presented can be extended to higher order. At second order one will need a population of M.s states 
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within a given cluster, and a population of A^c clusters of states. ITherefore the algorithm must follow a 
total population of A/'A^cA^s states. The problem will be to see if the resulting algorithm reaches accurate 
results within the numerically accessible values of AA, A^c, A^s- 

Some variants of the method should also be explored. In particular we have not exploited the variational 
formulation in the computation of the probability Q(h). One could also study in details the shape of the 
probability distributions of local fields in order to understand how they could parametrized in a simple but 
efficient way so that the free energy to be minimized does not involve a too large number of parameters. 

We thank Giulio Biroli, Cristiana Carrus, Enzo Marinari, Remi Monasson, Riccardo Zecchina and 
Francesco Zuliani for useful discussions and comments. 



APPENDIX I: WEIGHTED SUMS OF UNCORRELATED VARIABLES 



In this appendix we want to prove a useful little theorem, which applies to the computation of F^^' 
given in (^ and F^^) in (JH). 

Theorem: consider a set oi M{'> 1) iid random free energies /„, (a £ {1 . . . M}) distributed with the 
exponential density eq. (^ , and a set of M positive numbers a^. Then, neglecting terms which go to 
zero when M goes to infinity, the following relation holds: 



<'"(%^if^)'-<'°(?-)"4-(s?< 



(91) 



where {.)f denotes an average over the distribution of /. 

Corollary 1: In the same conditions as the theorem, for any set of M real numbers b^, one has: 



Corollary 2: In the same conditions as the theorem, for any set of M real numbers ba, one has: 



, J2aO'ablexp{-Pfa) 
T,a(^aexp{-Pfa) 



>/-( 



J2a(^abaexp{-/3fa) 

J2a O" exp(-/?/„) 



)f=X 



Y.c.<bl /Ea<^" 



EaflS 



EaOS 



(92) 



(93) 



Corollary 3: If the numbers Oq are M iid positive random variables, such that the average of a^ exists, 
which are uncorrelated with the fa, then one has: 



,, /^ Ea«aexp(-/3/a) Y _ /, (^^ V li r/x\N 

^ I E„exp(^/3/.) j ^f - ^1" 1^^ "'"""J ^f = -x ^" ^<""^'') 

where (.)a denotes an average over the distribution of a. 

Proof: we follow some of the techniques exposed in |^ . We start from the identity 

In I ^exp(-/3/a)aa J = / — exp(-i) - exp I -t^ exp(-/9/a)aQ J 

We choose to work with a regularised distribution of the M iid random variables /«: 

P{fa) = /3xexp(/3x(/„ - /,)) e{f, - f) , 



(94) 



(95) 



(96) 



where in the end we shall send M -^ cx), /c — > oo, with r = M cxp(— /?/c) fixed (the value of r is irrelevant). 
In this limit one has: 



(exp(-iexp(-/3/„)a„))/ ~ 1 - (to„)^ exp(-/3a;/c)r(l ~ x) , 
from which one deduces: 



(97) 
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(In I Y^ cxp(-/3/„)ao) \)f ^ — cxp(-i) - cxp(-r(l - a;)^ cxp(-/3a;/c) ^ 



^< 



X 



In ( r(l - xjf'M exp(-/3x/c) X! "S ) + -^— ^C* , (98) 



where C is Euler's constant. The quantity we need to compute involves subtracting the same expression 
with Oq substituted by one, which gives the desired result: 



Y^ exp(-/3/a)aa I - In I ^ < 

a. / \ a. 



(In ( l^ exp(~/3/„)a„ ]-^^{l^ cxp(-/3/„) \)f^i\nl^Y.<] (^9) 



The proof of corollary 1 is easily obtained by applying the theorem to the set of numbers aa exp(A&a), 
and taking the derivative with respect to A at A = 0. Similar generalised formulas can be obtained by 
taking higher order derivatives. The second derivative gives corollary 2. Corollary 3 is trivial. 

Remark: We notice that in the two limits a; — > and x ^ I the equations can be simply understood: 

• In the limit a; = 0, in a typical realization of the random free energies, only one weight w is equal to 
one and all the others are zero. Averaging over the realizations of free energies amounts to spanning 
uniformly the set of indices of this special non zero weight. 

• In the limit x — 1 the number of relevant w goes to infinity and each individual contribution goes 
to zero. An infinite number of term is present in the l.h.s. of eq. ( pi|) and the r.h.s. of the eq. ( pi| ) 
becomes In [(1/Af ) J2a '^"l' ^^ ^^ should. 
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